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Abstract 

We consider entropy solutions to the initial value problem associated with scalar nonlinear 
hyperbolic conservation laws posed on the two-dimensional sphere. We propose a finite volume 
scheme which relies on a web-like mesh made of segments of longitude and latitude lines. The 
structure of the mesh allows for a discrete version of a natural geometric compatibility condition, 
which arose earlier in the well-posedness theory established by Ben-Artzi and LeFloch. We 
study here several classes of flux vectors which define the conservation law under consideration. 
They are based on prescribing a suitable vector field in the Euclidean three-dimensional space 
and then suitably projecting it on the sphere's tangent plane; even when the flux vector in 
the ambient space is constant, the corresponding flux vector is a non-trivial vector field on 
the sphere. In particular, we construct here "equatorial periodic solutions", analogous to one- 
dimensional periodic solutions to one-dimensional conservation laws, as well as a wide variety 
of stationary (steady state) solutions. We also construct "confined solutions", which are time- 
dependent solutions supported in an arbitrarily specified subdomain of the sphere. Finally, 
representative numerical examples and test-cases are presented. 

Key words: hyperbolic conservation law, sphere, entropy solution, finite volume scheme, 
geometry-compatible flux. 
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1. Introduction 

In this paper, building on our earlier analysis in [6,2] we study in detail the class of 
scalar hyperbolic conservation laws posed on the two-dimensional unit sphere 

§ 2 = {(x,y,z) G M 3 , x 2 + y 2 + z 2 = 1} . 

We propose a Godunov-type finite volume scheme that satisfies certain important con- 
sistency and convergence properties. We then present a second-order extension based on 
the generalized Riemann problem (GRP) methodology [3]. 

It should be stated at the outset that an important motivation for this paper is the need 
to provide accurate numerical tools for the so-called shallow water system on the sphere. 
This system is widely used in geophysics as a model for global air flows on the rotating 
Earth [7]. In its mathematical classification it is a system of nonlinear hyperbolic PDE's 
posed on the sphere. Its physical nature dictates that it can be described "invariantly" , 
namely in a way which is independent of any particular coordinate system. Locally, it has 
the (mathematical) character of a two-dimensional isentropic compressible flow, whereas 
globally the spherical geometry plays a crucial role in shaping the nature of solutions 
-which, as expected for nonlinear hyperbolic equations, may contain propagating discon- 
tinuities such as shock fronts or contact curves. Thus, the relation of the present study to 
the shallow water system is analogous to the connection between Burgers' equation and 
the system of compressible fluid flow (say, in the plane). In fact, in light of this analogy 
it is somewhat surprising that in the existing literature so far, virtually all treatments, 
theoretical as well as numerical, were confined to the Cartesian setting. In particular, 
to the best of our knowledge, there have been no systematic numerical studies of scalar 
conservation laws on the sphere. 

Having introduced the scalar conservation law as a simple model for more complex 
physical systems, we should emphasize here also the intrinsic mathematical interest of 
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the model under consideration. It is already known (see [5] and references there) that even 
in the Cartesian setting, the two-dimensional scalar conservation law displays a wealth of 
wave interactions typical of the physical phenomena (such as triple points, sonic shocks, 
interplay of rarefactions and shocks coming from different directions and more). As we 
show here, "geometric effects", superposed on the (necessarily) two-dimensional frame- 
work, carry the scalar model still further. For example, the concept of "self-similar" solu- 
tions makes no sense here. In particular, one loses the Riemann solutions, a fundamental 
building block in many schemes (of the so-called "Godunov-type" ) . On the other hand, 
it allows for large classes of non-trivial steady states, periodic solutions, and solutions 
supported in specified subdomains. All these have natural consequences in developing 
numerical schemes; they offer us a variety of test-cases amenable to detailed analysis, to 
be compared with the computational results. 

In practical applications a finite volume scheme requires a specification of a coordinate 
system, where the symmetry-preserving latitude-longitude coordinates are the "natural 
coordinates" of preferred choice. The proposed finite volume scheme in this paper is based 
on these natural coordinates, but should pay attention to the artificial singularities at 
the poles. 

In [2] , a general convergence theorem was proved for a class of finite volume schemes 
for the computation of entropy solutions to conservation laws posed on a manifold. As a 
particular example, the case of the sphere S 2 was discussed, both from the points of view 
of an "invariant" formalism and that of an "embedded" coordinate-dependent formula- 
tion. In the present study we focus on the sphere § 2 and we actually construct, in a fully 
explicit and implemcntable way, a finite volume scheme which is geometrically natural 
and can be viewed as an extension of the basic Godunov scheme for one-dimensional 
conservation laws. Furthermore, we prove that our scheme fulfills all of the assumptions 
required in [2], which ensures its strong convergence toward the unique entropy solution 
to the initial value problem under consideration. We then describe the GRP extension 
of the scheme, whose convergence proof is still a challenging open problem. 

The theoretical background about the well-posedness theory for hyperbolic conserva- 
tion laws on manifolds was established recently by Ben-Artzi and LeFloch [6] together 
with collaborators [1,2,8]. An important condition arising in the theory is the "zero- 
divergence" or geometric-compatibility property of the flux vector; a basic requirement 
in our construction of a finite volume scheme is to formulate and ensure a suitable discrete 
version of this condition. 

We conclude this introduction with some notation and remarks connecting the present 
paper to the general finite volume framework presented in [2] . Following the terminology 
therein, we use an "embedded" approach to the spherical geometry, namely, we view 
the sphere as embedded in the three-dimensional Euclidean space K 3 . We denote by x 
a variable point on the sphere § 2 , which can be represented in terms of its longitude A 
and its latitude cf). Following the conventional notation in the geophysical literature we 
assume that 

< A < 2tt, --<<*<-, 

so that the "North pole" (resp. "South pole") is at <j> = f (resp. —(f) — |) and the 
equator is{</> = 0,0<A<27r}. (See Figure 1.) The coordinates in M 3 are denoted by 
(a;i, X2, Xs) G M 3 and the corresponding unit vectors are ii, i2, i3- Thus, at each point 
x = (A, <p) e S 2 , the unit tangent vectors (in the A, (f> directions) are given by 
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i\ = — sin A ii + cos A 12, 

i0 = — s'mcj) cos Aii — sin0 sinAi 2 +cos^>i3. 

It should be observed that while a choice of a coordinate system is necessary in practice, 
it always introduces singularities and the unit vectors given above are not well-defined 
at the poles and, therefore, in the neighborhood of these points it cannot be used for a 
representation of smooth vector fields (such as the flux vectors of our conservation laws). 
We also emphasize that the status of these two poles is equivalent to the one of any 
other pair of opposite points on the sphere. When such local coordinates are introduced, 
special care is needed to handle these points in practice, and this is precisely why we 
advocate a different approach. 

Continuing with the description of our "embedded" approach, we define the unit nor- 
mal n x , to § 2 at some point x by 

n x = cos (f> cos A ii + cos <f) sin A i 2 + sin 4> 13- 

Then, any tangent vector field F to § 2 is represented by 

F = F x i A + F - U 

and the tangential gradient operator is 

1 d d_ 

k cos <fi dX ' d<f> / 

Thus, the (tangential) gradient of a scalar function h(X, 4>) is given by 

^ , I dh . dh . .„ . 

COS (j) OA 0<j> 

and the divergence of a vector field F is 

V - F =cc^(^^ C0S ^ + ^)- (L2) 

Given now a vector field F = F(x, u) depending on a real parameter u, the associated 
hyperbolic conservation law under consideration is 

^ + V T -(F(x,«))=0, (x,()£§ 2 x[0,4 (1.3) 

where u = u(x, t) is a scalar unknown function, subject to the initial condition 

u(x,0)=uo(x), xe§ 2 (1.4) 

for some prescribed data u on the sphere. As mentioned above, we will impose on the 
vector field F(x, u) an additional "geometry compatibility" condition. 

An outline of this paper is as follows. In Section 2, we consider the construction of 
geometry-compatible flux vectors, while Section 3 is devoted to a description of several 
families of special solutions associated with the constructed flux vectors. In Section 4 
we discuss our (first-order) finite volume scheme, which can be regarded as a Godunov- 
type scheme. We prove that it satisfies all of the assumptions imposed on general finite 
volume schemes in [2], and we conclude that it converges to the exact (entropy) solution. 
In Section 5 we describe the (second-order) GRP extension of the scheme. Finally, in 
Section 6 we present a variety of numerical test cases. 
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Fig. 1. Web grid on a sphere 

2. Families of geometry-compatible flux vectors 

As pointed out in [2], every smooth vector field F(x, u) on § 2 can be represented in 
the form 

F(x, u) = n(x) x #(x, u), (2.1) 

where <fr(x, it) is a restriction to S 2 of a vector field (in M 3 ) defined in some neighborhood 
(i.e. , a "spherical shell" ) of § 2 and for all values of the parameter it. The basic requirement 
imposed now on the flux vector F(x,it) is the following divergence free or geometric 
compatibility condition: For any fixed value of the parameter 

V T -F(x,u)=0. (2.2) 

A flux vector F(x, it) satisfying (2.2) is called a geometry-compatible flux [6]. Note that 
this condition is equivalent, in terms of the nonlinear conservation law (1.3), to the 
following requirement: constant initial data are (trivial) solutions to the conservation 
law. In the case of the sphere S 2 the condition (2.2) can be recast in terms of a condition 
on the vector field <fr(x,it) appearing in (2.1). See [2, Proposition 3.3]. 

Our main aim in the present section is singling out two (quite general) families of 
geometry-compatible fluxes of particular interest, which arc amenable to detailed ana- 
lytical and numerical investigation. 
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The flux- vectors of interest are introduced by way of the following two claims. 

Claim 2.1 (Homogeneous flux vectors.) If the three-dimensional flux <&(x, u) = <J>(u) 

is independent of x (in a neighborhood ofE> 2 ), then the corresponding flux vector F(x,it) 
given by (2.1) is geometry-compatible. 

Proof. The following decomposition applies to any vector <!>(«) e R 3 in the form 

*(«) = h{u) h + f 2 (u) i 2 + h{u) i 3 , (2.3) 

so that F(x, u) = F\(\, <p, u) i\ + F^X, <fi, u) i^, with 

F x {\,<j),u)= /i(u) sin</> cos A + f 2 (u) sin</> sin A - f 3 (u) cos</>, 
F (A,(/),u) = -/i(m) sinA + / 2 (u) cos A. 

We can directly apply the divergence operator (1.2) to F(x, u) and the desired claim 
follows. □ 



Claim 2.2 (Gradient flux vectors.) Let h = h(x, u) be a smooth function of the 
variables x (in a neighborhood ofS 2 ) and u G R, and consider the associated three- 
dimensional flux 4>(x, u) = V/i(x, u) (restricted to x £ E> 2 ). Then, the flux vector F(x, u) 
given by (2.1) is geometry-compatible. 

Proof. We use the divergence theorem in an arbitrary domain D C § 2 with smooth 
boundary dD: 

/ V T • (F(x,u)) da = / F(x,v) ■ v(x)ds 

JD JdD 

(n(x) x V/i(x, w)) • i/(x) ds, 



idl 



IdD 

where i/(x) is the unit normal (at x) along dD C § 2 , da is the surface measure on § 2 , 
and (is is the arc length along dD. 

In particular, n(x) x v(x) = t(x) coincides with the (unit) tangent vector to 3D at x. 
It follows that the triple product (n(x) x V/i(x, u)) ■ v(x) = Vh(x,u) ■ t(x) is nothing 
but the directional derivative Van of h along dD. Since 



we thus find 



VoDhds = 0, 

dD 



/ V T • F(x, u) da = 0, 

JD 



and since this holds for any smooth domain D, we conclude that V T • F(x, v) = for all 

v e K. □ 



Remark 2.3 i. Claim 2.1 is a special case of Claim 2.2. Indeed, by taking in the latter 
h(x,u) — Xifi(u) + x 2 f2(u) + xzfz{u) we obtain the conclusion of the former. However, 
we chose to single out Claim 2.1 as a special case since it will serve in obtaining special 
solutions (Section 3) and in dealing with numerical examples (Section 6). 
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2. The steps in the construction of the gradient flux vector in Claim 2.2 are "linear 
in nature", namely if h(x,u) — hi(x,u) + h,2(x, u) then the corresponding (geometry- 
compatible) flux vectors satisfy F(x, u) = Fi(x, w) + F2(x, u). However, it is clear that the 
corresponding solutions to (1.3) do not add up linearly, due to the nonlinear dependence 
in u. 

3. Special solutions of interest 

3.1. Periodic equatorial solutions 

The scalar conservation laws discussed in this paper have two basic features: 

- The problem is necessarily two-dimensional (in spatial coordinates). 

- The geometry plays a significant role, inasmuch as the flux vectors are subject to 
geometric constraints. 

It should be noted that even within the framework of Euclidean two dimensional con- 
servation laws there is a great wealth of special solutions, displaying complex wave in- 
teractions, such as triple points, sonic shocks and more. We refer to [9,5] for detailed 
treatments of the theoretical and numerical aspects. 

In the situation under consideration in the present paper, geometric effects yield a 
large variety of non-trivial steady states, solutions supported in arbitrary subdomains, 
etc. In this section we consider such solutions by selecting some special flux vectors 
F(x, u) on § 2 . This is accomplished by making special choices of <fr(x,u) in the general 
representation (see (2.1)) F(x, u) = n(x) x <&(x, u), where <fr(x,u) is a restriction to S 2 
of a vector field (in M 3 ) defined in some neighborhood (i.e., "spherical shell") of S 2 and 
for all values of the parameter u. 

We begin our discussion with the case of periodic equatorial solutions, defined as 
follows. Taking f\{u) = f2{u) = in the general decomposition (2.3) so that, by (2.4), 

F x {X,^,u) = -f 3 (u) cos0, 
F 4> (X,4>,u)=0, 

the conservation law (1.3) takes the particularly simple form 

^-§x Mu) = ' (x,t)e$ 2 x{0,^). (3.1) 
In particular, obtain the following important conclusion. 

Corollary 3.1 (Solutions with one-dimensional structure.) Let u = u(X,t) be a 
solution to the following one- dimensional conservation law with periodic boundary con- 
dition 

Ou 

fo-d\fo&) = 0, 0<A<2tt, u(0,t)=u(27r,t), 

and let u = u(<f>) be an arbitrary function. Then, the function u(X, <f), t) = u(X, t) u(<j>) is 
a solution to the conservation law (3.1). 

It follows that all periodic solutions from the one-dimensional case can be recovered 
here as special cases. However, in numerical experiments the computational grid is two- 
dimensional, so it is not obvious that the accuracy achieved in the computation of the 
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former can indeed be achieved in the numerical scheme implemented on the sphere. This 
issue will be further discussed below, in Section 6. 

3.2. Steady states 

Let F = F(x, u) be a flux vector and uq : § 2 — > W. be an initial function such that 
V T • (F(x,u (x))) = 0. Then, clearly u is a stationary solution (or steady state) to the 
conservation law. In fact, we can show that there exist many (analytically computable) 
non-trivial steady state solutions, as follows. 

Claim 3.2 (A family of steady-state solutions.) Let h = h(x,u) be a smooth func- 
tion defined for all x in a neighborhood of§ 2 , and consider the associated gradient flux 
vector $ = V/i (as in Claim 2.2). Suppose the function uq : § 2 — > R satisfies the condition 

V y /i(y,u (x))| y=x = V x if(x), x e § 2 , (3.2) 

where H = H(x) be a smooth function defined in a neighborhood o/§ 2 . Then, u is a 
stationary solution to the conservation law (1.3). 

Proof. We follow the proof of Claim 2.2 and the notation therein. Using the divergence 
theorem in an arbitrary domain DCS 2 with smooth boundary dD, we obtain 

J V T • (F(x,m (x))) da = J F(x,u (x)) • u ds 

dD 

= j (n(x) x V x i?(x)) • i/(x) ds. 

dD 

where, as before, i/(x) is the unit normal, da the surface measure, and ds the arc length. 
In particular, n(x) x f(x) = t(x), the (unit) tangent vector to dD at x. It follows that 
the triple product (n(x) x V x i?(x)) • i/(x) = (V x _ff(x)) -t(x) is the directional derivative 
of H along dD. Thus, 

/ V T - (F(x,u o (x)))rfa = 0, 
Jd 

and since this holds for any smooth domain D, it follows that V T • (F(x,u (x))) = 0, 
which concludes the proof. □ 

The above claim yields readily a large family of non-trivial stationary solutions, as 
expressed in the following corollary. 

Corollary 3.3 Consider the flux vector F = F(x, u) given by 

F(x,w) = n(x) x (/i(u)ii), 

for an arbitrary choice of function f\ = fi(u). Then, any function u — u (xi) depending 
only on the first coordinate X\ is a stationary solution to the conservation law (associ- 
ated with this flux). In particular, in polar coordinates (A, <f>) any function of the form 
Uq(\,4>) — g(cos(/>cos A) is a stationary solution. 
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Proof. According to Claim 2.1 this flux vector is associated with the scalar function 
h(x,u) — Xifi(u). So we can invoke Claim 3.2 with H(x) = H(x\) such that H'(x\) = 
/i(uo(*i)). □ 

Remark 3.4 This corollary enables us to construct stationary solutions supported in 
"bands" on the sphere. This is accomplished by taking uq = uq(x\) to be supported in 
0<a<x\<P<l. Observe that this band is not parallel neither to the latitude curves 
(<fi = const) nor to the longitude curves (X = const). 

There is yet another possibility of obtaining stationary solutions, where all three coor- 
dinates are involved, as stated now. This example can also be derived from the previous 
one by applying a rotation in R 3 . 

Corollary 3.5 Consider the flux vector F = F(x, u) be given by 

F(x,u) = n(x) x ii + f 2 (u) i 2 + / 3 (u) i 3 ) 

= /(u)n(x) x (ii + i 2 + i 3 ), 

in which all three components coincide: fi(u) = faiu) = fo{u) = f(u). Then, any func- 
tion of the form Uo( x ) = Uo{%i + 22 + ^3), where uq depends on one real variable, only, 
is a stationary solution to the conservation law associated with the above flux. 

Proof. Following the proof of the previous corollary, we now take -ff(x) = H (xi+X2+x 3 ), 
where H' (£) = /(«o(0)- D 

Remark 3.6 In analogy with Remark 3.4, this result allows us to construct stationary 
solutions in a spherical "cap" (a piece of the sphere cut out by a plane). In Section 6 
below, we will provide numerical test cases for such stationary solutions. 

3.3. Confined solutions 

If in the conservation law (1.3) we have F(x, u) = for x in the exterior of some 
domain DCS 2 , identically in u G R, and if the initial function Uo(x) vanishes outside 
of D, then clearly the solutions satisfy u(x, t) = for x ^ D and alH > 0. We label such 
solutions as confined (to D) solutions. In view of equation (2.1) a sufficient condition for 
the vanishing of F(x, u) outside of D is obtained by 3>(x, u) = for x <^ D, identically in 
u E R. In view of Claim 2.2, this will follow if we choose h(x,u) such that h(x,u) ^ for 
x only in D. In particular, let tp = be a twice continuously differentiable function 
on R supported in the interval (a, /3) C (0, 1) and such that 3/3 2 > 1 and 3a 2 < 1. With 
an eye to computable test cases, we can use this function to generate solutions which are 
confined within the intersection of S 2 with the (three-dimensional) cube [a,/3] 3 . 

Claim 3.7 (A family of confined solutions.) Let ip be as above and let f = f(u) be 
any (smooth) function ofu^R. Define h = h(x,u) by 

h(x, u) = %l)(x 1 )ip(x 2 ) i>{x 3 ) f(u), 
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and let F(x, u) be the gradient flux vector determined in terms of h(x, u) as in Claim 2.2. 
Let DCS 2 be the spherical patch cut out from S 2 by the inequalities a < Xi < f3, i = 
1,2,3. Then, if the initial data Wo(x) is supported in D, the solution u = u(x,i) of the 
conservation law (1.3) associated with F(x, u) is supported in D for all t>0. 

Possible choices for a function ip : [a, 0\ —* R as in the claim are = sin 2 (fc£) for 
some integer k such that ka and kj3 are multiples of tt, or else tp(£) — (£ — a) 2 (^ — (3) 2 . 



4. Design of the scheme 

4.1. Computational grid 

The general structure of our grid is shown in Figure 1 , and its essential feature is the 
following. Every cell 1Z is bounded by sides which lie cither along a fixed latitude circle 
(4> = const.) or a fixed longitude circle (A = const.). We have 



K ■= {Ai < A < A 



2, 



<6< 



(4.1) 



as represented in Figure 2. In most cases, dlZ consists of the four sides of 1Z. However, 
across special latitude circles we reduce the number of cells, so that the situation (for a 
reduction by ratio of 2) is as in Figure 3. In this case the boundary dlZ consists of five 
sides, (so that the intermediate point (A3, 1^2) is regarded as an additional vertex), and 
even in this five-sided cell 1Z every side satisfies the above requirement. 



Phi ' 








^ Eprime D , 

Rectangle 




nu 


Phim - 






Lambdam 


Lambda 



Lambda 1 Lambda2 
Fig. 2. Rectangular cell TZ as part of grid on S 2 



The length of a side e = {Ai < A < A2, <fi = const.} equals (A2 — Ai)cos</>, while the 
length of a side e' — {4>i < <fi < (f>2, A = const.} is fa — <j>i- Consequently, the area A-ji 
of the cell 1Z is 



/■A2 rfc 

i n = dX 



cos(j)d(f) = (A 2 — A!)(sin0 2 — sin^i). 
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Fig. 3. Five-sided rectangular cell 1Z (on southern hemisphere of S 2 



4.2. Geometry-compatible discretization of the divergence operator 

Given any rectangular domain 1Z of the form (4.1), the approximate flux divergence is 
now derived as an approximation of the integral of the flux along the boundary dlZ, 
divided by its area, as follows: 



where ds is the arc length along dlZ and u is the outward-pointing unit normal to 
dlZ C § 2 . In the limit A 2 ,^2 —> Ai,^i the approximation (4.2) to the divergence term 
approaches the exact value (1.2). 

We need to check that the geometric compatibility condition (2.2) is satisfied for the 
approximate flux divergence. This requirement will be taken into account in formulating 
our finite volume scheme for (1.3). 

Consider now the actual evaluation of the term In defined in (4.2) and consider the 
cell shown in Figure 2, under the assumption that u = u(X, (p, t) is smooth on 1Z. We 
propose to approximate the flux integral along each edge of 1Z in the following way. As 
in Section 2, let us decompose the flux into its (A, 0) components: 



On each side the integration is carried out by (i) taking midpoint values of the appro- 
priate flux component, and (ii) using the correct arc- length of the side. We designate 
the midpoints of the edge e as A e ' m = (Ai + A 2 )/2 and <p e ' m — <\>\ (see Figure 2), and 
likewise for the edge e'. 

Throughout the rest of this section we restrict attention to the gradient flux vector 
constructed in Claim 2.2. In particular, it comprises the class of homogeneous dux vectors, 




(4.2) 



dTZ 



F(x, u) = -F\ (A, <j>, u)i\ + i^(A, </>, u)i . 
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given by (2.3)-(2.4). 

Taking u as constant u — u e - m along the side e € dlZ, the total approximate flux is 
given by 

-i approx 

F{x,u)-uds =-(h{e 2 ,u e ' m )-h{e 1 ,u e ' m )), (4.3) 

where e x ,e 2 are, respectively, the initial and final endpoints of e (with respect to the 

sense of the integration). 

Summing up over all edges we obtain: 

Claim 4.1 (Discrete geometry-compatibility condition.) Consider the gradient flux 
vector constructed in Claim 2.2. Then, if u = const., 1-r = 0, so that 

[V T -F(x,u)]° ppr ° x = 0, 

and thus a discrete version of the divergence- free condition (2.2) holds. 

Remark 4.2 The claim above applies to gradient flux vectors in Claim 2.2, and, in par- 
ticular, to homogeneous flux (2.3) -(2.4). On the other hand, for a more general geometry- 
compatible flux F(x, u), such a result can be obtained only if the dependence on x is 
integrated exactly along each side, a requirement that must be imposed on the scheme. 

4.3. Godunov-type approach to the numerical flux 

We continue to deal with the gradient flux given in Claim 2.2. We assume different 
(constant) values of u = u(A, <f>, t) in grid cells and evaluate the numerical flux values at 
each edge from the solution to a Riemann problem with data comprising these values 
u(X, (j), t) in the cells on either side of that edge. At the midpoint (A e,m , e,m ) of each side e 
we solve the Riemann problem in a direction perpendicular to e, and denote the resulting 
solution u e - m . The corresponding fluxes are then evaluated as F(A e ' m , (j) e ' m ,u e ' m ). 

We can split Eq. (1.3) by invoking the explicit form of the divergence (1.2), getting 
du 1 d 

+ 7TJT F a(A,0,w) = O for the side e' : A = A 2 , (4.4) A 

at cos a\ 

ttt —7 ^i( F 4>{\ 4>, u)cos0 N ) = for the side e : <f> = </> 2 , (4.4),* 

at cos <po<p\ / 

Consider two adjacent cells, as in Figure 4 or in Figure 5. By fixing <f> = <fi e ' m (resp. A = 
A e ' m ) in (4.4)a (resp. (4.4)^ ) we can evaluate u 
\ = \ e m (resp. 4> = 4> e > m ). 



We include here some remarks that will be useful in the implementation of the scheme. 
Consider an homogeneous flux vector as in Claim 2.1 so that its components are given 
by (2.4). Suppose that u(X, <fi, t n ) — u L (resp. u(\,cj>,t n ) = u R ) in the cell {Ai < A < 

A 2 , 0i < 4> < 4>2} (resp. {A 2 < A < A3, <f>i < 4> < 2 j), as in Figure 4. At the point 

Af(A e '' m ,0 e '' m ) Eq. (4.4) A takes the form 

^ + tan ¥' m ^~ X cos A + sin A ) - = °- ( 4 - 5 ) 
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Lambdal Lambda2 Lambda3 

Fig. 4. Two A-adjacent cells with constant states u L , ur 
Phi3 



Phi2 




u=u_L 



Phil I | I 

Lambdal Lambdam Lambda2 

Fig. 5. Two (^-adjacent cells with constant states u L , Ur 

Setting 

ff (A,u)=tan^ m (/i(u)cosA + / 2 (u)sinA) - f 3 (u), (4.6) 

we see that equation (4.5) is the scalar one-dimensional conservation law 

^ + ±g(X,u) = 0, t>t n (4.7) 

subject to the initial data u = u L (rcsp. u = u R ) for A < A 2 (rcsp. A > A 2 ). 

Likewise, we repeat the former analysis for ^-adjacent cells by taking the constant 
states u(X, <f), t n )=u L , u(X, <f>, t n )—u R in cells {Ai<A<A 2 , 4>i<4> < ^2}, {Ai<A<A 2 , <f)2<4><4> 
as depicted in Figure 5. At the point M(A = X e ' m , (f> = <f> 2 ), the equation (4.4)0 then takes 
the form 

^ + ^-|rf-sinA e ' m cos0/ 1 ( U ) + cosA e ^cos ( /)/ 2 ( M )) =0. (4.8) 
01 cos <po<p\ / 
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We then set the 0-flux function 

k(<j>,u)= (-sinA e ' m /i(u) + cosA e < m / 2 (u))cos0 (4.9) 

so that equation (4.8) is the scalar one-dimensional conservation law 
du ' d 
at cos (f> oq> 

subject to the initial data u — u L (resp. u — u R ) for (f> < fa (resp. <f> > fa). 

4.4. Solution to the Riemann problem 

The solution at the discontinuity A=A 2 at the initial time t = t n is given by the 
Riemann solution to (4.4) \ . For simplicity of the presentation we specialize here to the 
flux (4.7). Since the dependence of g(X,u) on A is smooth, this solution is obtained by 
fixing A = A 2 , thus solving the classical conservation law 

^ + ^.g(A 2 ,«) = o, t>t n (4.11) 

subject to the initial jump discontinuity of u. 

We denote this solution by u 2,m . Observe that the flux g(X,u) in (4.11) is in gen- 
eral non-convex. The Riemann solution may therefore consist of several waves. It is a 
self-similar solution depending only on the slope (A — A 2 )/(i — t n ). The value u 2 ' m is 
the value along the line A=A 2 . It therefore corresponds cither to a sonic wave, namely 
</(A 2 , u 2,m )=0, or to an "upwind value" u—u L (resp. u—u R ) in the case where all waves 
propagate to the right (resp. left). 

Actually, the procedure for solving the Riemann problem in the case of a nonconvex 
flux function g(A 2 ,u) is well-known and goes back to classical works by Oleinik and 
others. We recall it here briefly. Assume first that u L <u R . Consider the convex envelope 
of g, namely, the largest convex continuous function g c , over the interval [u^uj, such 
that g c <g at all points. Clearly, g c =g in "convex sections" of the graph of g, while it 
consists of linear segments when g c <g. It is easy to see that the "convex segments", 
where g—g c , represent rarefaction waves (in the full Riemann solution) while the linear 
segments represent jumps (i.e., shock waves). In particular, the solution u 2 ' m is given by 
the following formula: 

where g(X 2 , v min ) < g(X 2 , v) for all v&[u L ,u R }. (4.12) 

There are in fact three possibilities for this solution: 

a) u L < u 2 ' m < u R , which implies that </(A 2 , u 2 ' m ) — (a sonic point). 

b) u 2 ' m = u L , the whole wave pattern moves to the right. 

c) u 2 ' m = u R , the whole wave pattern moves to the left. 



Similarly, in the case u L >u R , we construct the "concave envelope" of g, namely, the 
smallest concave continuous function g c such that g c >g. Again the linear segments cor- 
respond to jump discontinuities while the concave segments (g=g c ) correspond to rar- 
efaction waves. The solution to the Riemann problem is now given by u 2 ' m ~v max , where 
ff(A 2 , v max )>g(X2, v), v € [u R ,u L ]. As above, there are three possibilities for the solution 
(sonic, left-upwind, or right-upwind). 
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Replacing in the foregoing analysis the A- flux function g(\ 2 ,u) by the fflux function 
k((f>2,u), the equation (4.10) reads 

^ + ^(-sinA^AH +cosA 2 ^/ 2 (w)) =0, t > t n . (4.13) 

We get the Riemann solution to (4.13) in the three cases a), b), c) as above. 

4.5. Convergence proof 

The computational elements ("grid cells") are denoted in [2] by K. Their sides are 
denoted by e and the flux function across e is given by / e> jf (w, v), where u is the (constant) 
value in K and v is the value in the neighboring cell (sharing the same side e) K e . In 
our grid of the sphere, some cells arc actually pentagons; these arc the cells whose lower- 
latitude side (along a latitude (f> = const) borders the two higher-latitude sides of the two 
lower-latitude neighbor cells, as shown in Figure 3 for the southern hemisphere grid. For 
such cells, the lower-latitude side consists of two faces, each one of them common with 
one of the lower- latitude neighboring cells. 

With this construction of the grid, we can check the conditions in [2] imposed on the 
numerical flux. It is important to keep in mind that we are dealing with the gradient flux 
vectors given by Claim 2.2. 

Claim 4.3 (Convergence of the proposed scheme.) Consider the first- order finite 
volume scheme described above. Assume that the flux vector has the gradient form in 
Claim 2.2. Let f e .n{ u , v ) be the numerical flux calculated on the side e of the computa- 
tional celllZ, using (4.3), where the midpoint value of u is obtained from the Riemann 
solution. Then f e ^i(u, v) satisfies the assumptions (5. 5)- (5. 7) of [2], and the numerical 
solution converges to the exact solution as the maximal size of the grid cells shrinks to 
zero. 

Proof. Consider the flux across a longitude side e : A = A 2 , which is given by F\ in 
the equation (4.4)a . The procedure for integrating the flux across e is described by 
(4.3), while in Subsection 4.4 the calculation of F\(A 2 , </> 2 ' m , it 2 '" 1 ) is described. It can be 
summarized as follows. 

First, the solution u 2 ' m to the Riemann problem associated with equation (4.4)a is 
found, assuming u,v to be the values on the two sides. However, note that F\ depends 
explicitly on <j), and to be precise we need to replace in (4.4) a the mean value 2,m by (f>. 
Thus, we find u 2 ' m = u 2 > m (<f>). 

Clearly, in the case M = nwe get identically u 2 ' m {<p) = u = v and so the exact flux 
satisfies 

F x = F x (\ 2 ^,u 2 > m ) 
and its integration will give exactly the approximate value 

/»,*(«, v) = -(h(e 2 ,u 2 ' m ) - h(e\u 2 ' m )), 

as in (4.3). Thus, condition (5.5) in [2] is satisfied. 

Clearly, the conservation property (5.6) is satisfied even with the approximate defini- 
tion. 
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Also, the flux as defined in (4.3) makes it easy to check (5.7), as the flux is independent 
of <f> and the monotonicity is thus a result of general properties of the Riemann solver 
(even for nonconvex fluxes). For example, if u < v, one considers the convex envelope of 
F\, as defined in (4.4) a (with <j> = <f> 2 ' m ) and then considers u 2 ' m as the minimal value 
on this envelope (over [u, «]). Clearly changing u upward will cither change u 2 ' m upward 
or leave it unchanged. This completes the proof. □ 

5. Second-order extension based on the GRP solver 

To improve the order of accuracy, we consider again the cell Ai<A<A 2 , </>i<0<0 2 and 
assume that u is linearly distributed there. We use u L ^\, u L ^ (resp. u Rj \, u R ^) to denote the 
slopes in the cell to the left (resp. right) of the side A=A 2 . We also denote by u L ((p) (resp. 
u R (<fi)) the limiting value (linearly distributed) of u at A=A 2 — (resp. A=A 2 +). Clearly, 
the solution to the Riemann problem across the discontinuity is a function of <f>, and 
we denote it by u 2 ' m ((j)), which conforms to our notation in Subsecion 4.4 above (where 
u was constant on either side of the discontinuity). The value of u 2 ' m (<p) is obtained 
by solving the Riemann problem associated with Eq. (4.4)> with (f) 2 > m replaced by <fi, 
subject to the initial data u L (<p), u R ((f>). Restricting to the middle point </> = <f) 2 ' m , the 
solution u 2 ' m ((j) 2 ' m ) (at A = A 2 '" 1 ) is in one of the three categories listed above (i.e., sonic, 
left-upwind, right-upwind). By continuity, the solution u 2 - m (<f)) will still be in the same 
category for — (f> 2 ' m sufficiently small. The solution at (A 2,m ,0 2 ' m ) varies in time and 
the GRP method deals with the determination of its time-derivative at that point. 

Accounting for the variation of the solution over a time interval enables us to modify 
the Godunov approach to the determination of edge fluxes , as presented in Section 4.3. 
We assume that the flux vector depends explicitly on x, as in (2.1). In what follows we 
use for simplicity the "imbedded" notation x = (xi, x 2 , x 3 ) for a point on the sphere (see 
the Introduction), along with the corresponding spherical coordinates A, <p. We further 
assume that the vector field 3> is given by the following extension of (2.3) 



The zero-divergence identity is obtained as a result of expressing $ as a gradient Vh in 
the sense of Claim 2.2. 

For our choice of <I> such a representation of <I> as gradient of h is obtained when h is 
taken as 



and qj(xj) = r'^xj), j = 1,2,3. 

Using (1.2) together with the geometry-compatibility property, we get an explicit form 
of the conservation law (1.3) in our case as 

^" -sm\q 1 (x 1 )-^-f 1 (u) + cos\q 2 (x 2 )-^-f 2 (u) 



The numerical approximation to this equation requires an operator splitting approach, 
where the derivatives with respect to (f> and A are considered separately. We note that 



$(x, u) = V x /i(x, u) 

= ii + q2(x 2 )f2(u) i 2 + q3(x 3 )f 3 (u) i 3 



(5.1) 



h(x,u) = ri(a;i)/i(u) + r 2 (x 2 )f 2 (u) + r 3 (x 3 )f 3 (u) , 



(5.2) 




(5.3) 



1G 



such a splitting has already been implemented in the Godunov case, (4.4), in the most 
general case. In that case, no use has been made of the geometry-compatibility property. 
Indeed, this has no bearing on the first-order scheme since the solution to the Riemann 
problem is obtained by "freezing" the explicit dependence on A, <fi (and, in particular, 
ignoring the terms involving the derivatives with respect to this explicit dependence). 
In the present (second-order) situation we proceed as follows. 
The "A-split" equation obtained from (5.3), is 

— + tan^ 2 ' m (^i(2:i)cosA — fi(u) + q 2 (x 2 ) sin A — ,f 2 (u)j - q 3 (x 3 ) g^f3( u ) = °- 

(5.4) 

Note that the coefficients are retained as functions of A and are not "frozen" at A = 
A 2,m . This is of course due to the fact that in employing the GRP scheme we consider 
A-derivatives on either side of the edge, so as in any limiting analysis, we must first let 
A -> A 2 < m , then substitute A = A 2 '" 1 . 

The A-edge flux function g(X,u) (compare (4.6)), is now extended to g(x,u) as 

g(x,u) = tan0 2 ' m ^(7i(a;i) cos A/i(u) + q 2 (x 2 ) sin A/ 2 (u)) - q 3 {x 3 )f 3 (u), 

and the scalar one-dimensional conservation law under consideration is now rewritten as 
an equation with a source term (a balance law) 

du d , . 

~dt + <9A' 9 ^ X ' U ' = ' 1 > *" 

S\ = tan^ 2 ' m (^/i(M) — (qi(xi)cosA) + f 2 (u) — (q 2 {x 2 )sin\)j - f 3 (u)—q 3 (x 3 ), 

(5.5) 

subject to the initial data (for u and its slope) u L ((p 2 ' m ), u L ^\ (resp. u R ((f> 2 ' m ), u Ri \) for 
A < A 2 (resp. A > A 2 ). Observe that the equation is written in a "quasi-conservative 
form" , which offers more convenience in the GRP treatment [3, Chap. 5]. The right-hand 
side term S\ is just the result of the A differentiation of the flux g(x, u). Obviously, 
the geometry-compatibility condition implies that this source term should cancel out 
with the corresponding source term in the "0-split" equation. The solution u 2,m to the 
Riemann problem is obtained by freezing the coordinate A at its edge value, so that, in 
particular, the source term in (5.5) can be taken as zero at this stage. 

In the framework of the GRP analysis, the source term S\ is added to terms arising 
from the piecewise-linear initial data, in producing the time-derivative of the solution 
u 2 ' m {4> 2 ' m ) + ^(\ 2 > m ,(t> 2 > m ,t n +)^f, At = t n+1 - t n . As explained above, u 2 ^ m (^ m ), 
the solution to the associated Riemann problem, is obtained by using the "edge val- 
ues" Ul(4>), u r (</>). It remains, therefore, to determine the instantaneous time-derivative 
|f (A 2 ' m ,<?!) 2 ' m ,i„+), as is outlined below. 

The time-derivative of u is given by 

^(\ 2 ' m ,(l> 2 ' m ,tn + ) = -U m .X J^ff(x,u)| A 2,m i0 2,m iU 2,m, 

where the slope value w TOj a is obtained by "upwinding", determined by the associated 
Riemann problem as follows (we start with the "easy" categories b) , c) above) . 
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b) u 2 ' m = u L ((f) 2 ' m ). Then, the wave moves to the right and we set 

u m ,\ = u L ,\. 

c) u 2 ' m — u R {<p 2 < m ). Then, the wave moves to the left and we set 

Urn, A = U R ,\. 

Finally, the first category deals with the sonic case. As noted above, it remains sonic 
in the neighborhood of 2,m , so that we have there J^g(x, u) \\2,™ ^2,™ jU 2, m . The time- 
derivative of u reduces therefore to 

|^(A 2 ,0 2 <"\i=t„+) = O. 
at 

Finally, the "0-split" equation obtained from (5.3), is treated in analogy with the 
"A-split" procedure outlined above. 

6. Numerical tests 

6.1. First test case: equatorial periodic solutions 

Here, the conservation law takes the form (3.1) and the flux function and initial data 
are given by 

fi(u) = h{u) = 0, h{u) = -2tt (u 2 /2), 

'sin A, < A < 2tt, < (f)< tt/12, (6.1) 

0, otherwise. 



u(A,0,O) 



As discussed in Section 3 (see the discussion of solutions to (3.1)) it is clear that the 
solution here (as a function of A) is identical to the periodic solution for the Burgers 
equation in Mr, with periodic boundary conditions on [0 < x < 2ir\. However, we compute 
the numerical solution here on our spherical grid, and we need to check not only that 
it conforms with the one-dimensional case but that it does not "leak" beyond the band 
supporting the initial data. The results at the shock formation time t s = l/2n are shown 
in Figure 6 for AA = 2tt/16, in Figure 7 for AA = 2tt/32 and in Figure 8 for AA = 2tt/64. 
These GRP solutions to (4.7) clearly converge to the exact solution with refinement of 
the A grid, and are comparable to the corresponding solution to the scalar conservation 
law in R 1 with Ax = 2?r/22. 



6.2. Second test case: steady state solutions 

We refer to Corollary 3.3 and using the notation there we take the flux vector and 
initial data as: 

h(u)=u 2 /2, / 2 («) = / 3 («) = 0, 
u(X, 4>, 0) = cos A cos 4>. 

Using the terminology of Corollary 3.3 we sec that the initial function is the "simplest" 
possible function, corresponding to g(x\) = x\. 
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0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 



Fig. 6. Exact, GRP/SCL and GRP/SPHERE (AA=2tt/16) solutions to the IVP (6.1) at t = l/2n 




0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 



Fig. 7. Exact, GRP/SCL and GRP/SPHERE (AA=2tt/32) solutions to the IVP (6.1) at t = l/2n 



As is shown in Figure 9, the numerical solution remains nearly unchanged in time 
after being subjected to integration up to t = 5 by the GRP scheme with constant time 
step At = 0.05, the color maps of u(X 7 (j) 7 t) at the initial and final times are virtually 
indistinguishable. The shown grid has latitude step A0 = n/60, and an equatorial 
longitude step AA = 7r/128. A measure Udiff to the numerical solution error is defined 
as the area-weighted difference |u(A, <f), 5) — u(X, </>,0)|, obtained by summation over all 
grid cells. In this case we obtained u^iff — 0.0093, which is small relative to the full range 
Umax — u m i n = 2. Hence, the GRP scheme produces an approximation to the steady- 
state solution u(A, 0, t) — u(X, 4>, 0) over § 2 . This test case demonstrates that the scheme 
computes correctly the time-evolution for the non-constant data (6.2), by calculating an 
approximately zero value for the flux divergence in computational cells. 
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0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 
Fig. 8. Exact, GRP/SCL and GRP/SPHERE (AA=2tt/64) solutions to the IVP (6.1) at t = 1/2tt 




Fig. 9. Steady-state initial data (and solution) to the IVP (6.2) at t = 5. 
Color map range scaled to {^mini'^ma.x ) = (-0.998,0.998). 

6.3. Third test case: confined solutions 

We take (as in Claim 2.2) 3?(x, u) = V/i(x, it), where /i(x, u) = ip(xi)xifi(u). The 
function ip{ x i) is defined by 
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1, x x < 0, 



0. 



6x1 - 
2 



o 3 

h 71 a;i ' 



0<xi< 



(6.3) 



The flux vector is then given by 

F(x, u) — n(x) x 3>(x, u). 

The solution is clearly confined to the sector x\ < ^ of the sphere. Its boundary is a 
circle which intersects the meridian A = at cf> — ~ . 
The flux in the subdomain x\ < is given by 

F(x,u) = n(x) x 

so if we take the initial data as ijj(xi)uo{x{) , where uq is the steady state solution of 
the second test case (and also the same fi(u)), the solution remains steady in that part, 
namely, in x\ < 0. Clearly, it evolves in time in the region < x\ < but vanishes 
identically (for all time) if ^ < x\. 

The confined IVP was integrated in time up to t — 5 by the GRP scheme, using the 
same grid and time step as in the second test case (Subsection 6.2). The solution is 
represented by the color map in Figure 10. Comparing it to the corresponding initial 
map (not shown here), it seems nearly unchanged. In fact, the initial-to-final difference 
measure obtained is Udiff = 0.0057, which indicates a nearly steady solution in the strip 
< xi < y/l/2 . This test case demonstrates that the scheme computes correctly the 
time-evolution for the non-constant "confined" data (6.3). 
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